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Abstract 

In this work, we systematically generalize the Evans function methodology to address vector 
systems of discrete equations. We physically motivate and mathematically use as our case example a 
vector form of the discrete nonlinear Schrodinger equation with both nonlinear and linear couplings 
between the components. The Evans function allows us to qualitatively predict the stability of the 
nonlinear waves under the relevant perturbations and to quantitatively examine the dependence of 
the corresponding point spectrum eigenvalues on the system parameters. These analytical predictions 
are subsequently corroborated by numerical computations. 
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1 Introduction 



In the last few years, the study of solitary waves in multi-component systems has drawn a large focus of 
attention both in the physics of Bose-Einstein condensates [1, 2, 3, 4, 5, 6, 7], and in that of nonlinear 
optical fiber and waveguide arrays [8, 9]. 

In the case of BEC dynamics, if the condensates are in a deep optical lattice, then their dynamics 
can be described by a discrete nonlinear Schrodinger (DNLS) equation [10, 11, 12]. Multi-species 
condensates can arise, in this setting, as mixtures of different spin states in 87 Rb [13, 14] and 23 Na 
[15] condensates. Furthermore, theoretical studies have discussed the possibilities of mixtures between 
different atomic species, such as Na-Rb [16, 17], K-Rb [18, 19], Cs-Rb [20] and Li-Rb [21]. A two- 
species BEC, in fact in a 41 K- 87 Rb mixture, has been recently reported [22]; finally, a mixture of 
7 Li and 133 Cs was also experimentally investigated [23] (albeit in a non-condensed state). In the above 
contexts, the relevant theoretical models, in the presence of an optical lattice trapping [4], consists of one 
DNLS equation per atomic species. These equations are coupled by nonlinear cross-phase modulation 
(XPM); however, linear coupling terms can also be applied to the description of BEC dynamics in this 
mean-field approximation. The nonlinear interaction between the components is generated by (inter- 
species) atomic collisions, while linear coupling may be readily induced by an external microwave or 
radio-frequency field which induces Rabi [24] or Josephson [25] oscillations between populations of the 
two states. 

Numerous issues have been considered, especially in the continuum (i.e., non-lattice) counterpart of 
the above framework. Among them, are ground-state solutions [16, 26, 27], small-amplitude excitations 
[17, 28], formation of domain walls (DWs) between immiscible species [29, 30, 31, 32], bound states 
of dark-bright [33] and dark-dark [34], dark-gray, bright-gray, bright-antidark and dark-antidark [35] 
complexes of solitary waves, spatially periodic states [36] and modulated amplitude waves [37] among 
others. 

A realization of the above discussed system is quite possible also in the setting of the coupled 
optical waveguides. In that case, the two species represent either two orthogonal polarizations of light, 
or two signals with different carrier wavelengths. In the latter case, the linear coupling can also be 
implemented, by twisting the waveguides or elliptically deforming them, in the case of linear or circular 
polarizations, respectively. In this context, strongly localized vector (two-component) discrete solitons 
have been identified, both in one-dimension [38, 39], as well as in two dimensions [40]. Similarly to 
their continuous counterparts [41], these vector solitons may have components of different types (bright, 
dark, or antidark). In particular, symbiotic bright-dark and dark-antidark pairs were predicted in such 
systems [38, 39]. Such solitary waves have also been observed experimentally, see e.g., [42]. 

These recent developments, on both theoretical and also experimental aspects of such multi-component 
systems underscore the relevance of further systematic studies of their coherent structures and, impor- 
tantly, also of their stability. More specifically, in the discrete setting, the stability of solitary waves 
in (a single-component) DNLS and related models was addressed in [43]. However, to the best of our 
knowledge, there were no works developing techniques systematically addressing the stability of waves 
in the context of vector discrete equations. This formulates the problem that the present study aims 
at addressing. In particular, we use as our starting point a vector integrable discretization of [44, 45] 
(see also [46] and references therein). We then include perturbations, through the inclusion of physi- 
cally relevant terms, such as ones that appear in the framework of the coupled DNLS models or the 
linear coupling discussed above (see e.g. [24, 25]). We focus on the influence of such terms on the 
point spectrum eigenvalues of the linearization around the perturbed solitary waves. To monitor the 
latter, we generalize the Evans function methodology of [43] (see also [47] for the continuum setting) 
to address the vector case and obtain good agreement of the relevant predictions with the qualitative 
phenomenology and the quantitative dependence on (the perturbed) system parameters. 
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Our presentation will be structured as follows. In section 2, we are going to give the general setup 
and notation of the coupled equation system. In section 3, we are going to develop the Evans function 
for the vector case and use it to obtain the eigenvalues in the case of DNLS-like perturbations. In 
section 4, we are going to give the corresponding results in the presence of linear coupling. In these two 
sections, the analytical results will be complemented with numerical computations. Finally, in section 
5, we summarize our findings and present our conclusions. 

2 The CDNLS equation and Notation 

We consider the case of two coupled discrete NLS (CDNLS) equations in the form 

iu n = -A 2 u n - (\u n \ 2 + \v n \ 2 )(u n+1 + u n -i) - 2eu n (\u n \ 2 + \v n \ 2 ) 

iv n = -A 2 v n - (\v n \ 2 + \u n \ 2 )(v n+ i + v n -i) - 2ev n (\v n \ 2 + \u n \ 2 ) (2.1) 

with 

A 2 := j^{u n +i ~ 2u n + u n -i) 

; h is the step-size of discretization and n is the lattice site index. Our motivating physical setting stems 
from the consideration of coupled hyperfine states of an atomic species (such as the spin states |1, — 1 > 
and |2, 1 > of 87 Rb; see [13, 14] for relevant details), examined in the presence of a quasi-one-dimensional 
deep optical lattice potential. A vector generalization of the derivation of [12] would establish the 
coupled discrete NLS model (with near- unity nonlinearity coefficients [14]) as the appropriate underlying 
mathematical framework. However, a starting point that is more tractable analytically is the e = 
variant of the above model of Eqs. (2.1) (i.e., the "Ablowitz-Ladik limit" [44, 45]). We therefore study 
first the existence and stability of discrete solitons for CDNLS equations starting from the e = limit 
and using continuation in e small to examine the persistence of such solutions in the non-integrable limit 
with the onsite nonlinearity. We study analytically the stability of discrete solitons of the above coupled 
system for C = 1/h 2 = 1 through the Evans method and the properties of reduced discrete dynamical 
system (perturbed 4-dimensional mapping). Based on the integrability of this map one can provide 
an analytic expression of the discrete Evans function. Concerning small perturbations of integrable 
CDNLS model, we wish to study the stability criteria for discrete solitons and test these predictions 
against numerical simulations. We note in passing that this type of discretization does have the relevant 
Manakov continuum limit as h — > 0, where the equations become 

iut = -An- (1 + e) (ju| 2 + |v| 2 ) u, (2.2) 
iv t = -Av- (1 + e) (|u| 2 + H 2 ) v, (2.3) 

up to rescaling of the amplitudes (u — ► u\J\ + e and v — ► vy/1 + e) or one of space and time (i.e., 
x — > xy/1 + e and t — > t(l + e)). Notice that this particular discretization is motivated by the wide 
volume of work on the so-called Salerno model [48] which is also considering a mixed nonlinearity of 
combined local and nearest-neighbor nonlinear term. 

We are interested in stationary solutions of the system of equations (2.1) and make the ansatz 

Un{t) = q n e~ iuJt , v n (t) = Pn e- iujt , (2.4) 

with a uniform rotation frequency to. Using Equations (2.4) in Equations (2.1), we arrive at two 
coupled second-order difference equations 
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Qn+1 + Qn-1 
Pn+1 +Pn-1 



2 "^ 2 q _ £2h 2 q {ql+Pl) 
2/2 2^ " n yr! 



For e = 0, we obtain the integrable standard-like map with two invariants: 



h = h 2 



h = 



( <7n + qI-1 - 2c QnQn-l ) + (pi + Pn-1 ~ ^PnPn-1 ) + h 2 (<7n + Pl){Qn-l + Pn-1 



h 2 



Qn + Qn-1 ~ 2c <7n<7n-l + P n + Pn-1 ~ ^PnPn-\ H ( ^n-l + PnPn-1 J 



where c = 2 — oj/i 2 . 
Upon setting 



= (Qu ~ Qn-l)/h, S n = (p n - p n -l)/h, 

the steady-state problem for CDNLS can be formulated as 



(2.5) 



Qn+l 


= (cx - 


l)?n + r n h + e( 


-2h 2 )q n {ql+p 2 n )x 


r n +i 


= (cx - 


l)|+r n + £ (- 


- 2h)q n (q 2 n + P 2 n )x 


Pn+l 


= (cx - 


l)Pn + Sn/l + £( 


-2h 2 ) Pn (ql+p 2 n )x 


s n+l 


= (cx - 


i \ £* n , / 
l)y + S n + £( " 


- 2h )p n (q n + Pl)X, 



where 



X 



(2.6) 



The equations (2.6) are written so that when e = 0, they are then exactly the steady-state problem for 
integrable CDNLS. 

When e = the solutions are given by 



Qn(0 
Pn(0 



aisinh2H / 

\A*i + a 2 
a2sinh2VU 



sech(TFn + £), 



sech(VFn + £), 



(2.7) 



with W = cosh _1 (c/2) = 2W. For small h, W « y/—ujh. Note that these solutions precisely describe 
the stable and unstable manifolds of the fixed point (0,0,0,0), and that these manifolds intersect non- 
transversely. This is a non-generic phenomenon for standard-like maps, and hence it is expected that 
the intersection, if it persists, will be transverse for e > 0. 

Addressing the linear stability of the wave, we set u n = u R + iUj , v n = v R + Wj and linearizing 
CDNLS about the wave, one has the linearized problem 

d t u R = -L_u I 

dtUj = L Q+ u R + L Q v R 



d t v R = -L_v z 

d t Vj = L p u R +L p+ v R 

Hence, the operators L_, L Q1 L p , L Q+ , L p+ are given by 



L_ 


= A 2 + co + (Q 2 n + PZ)(e d 


+ e" 9 ) 


L Q+ 


= A 2 + co + 2{Q n+1 + Q n . 


-i)Q„ + (Q f 2 t + P n 2 )(e 9 + e- 9 ) 




= 2Pn(Qn+l + Qn-l) 




L p+ 


= A 2 + u + 2{P n+1 + P„_ 




L p 


= 2Q n (P n+ i + P n _i) 





where e ±d u n = u n ±\. Upon setting 



it is not difficult to check that 

/ du>Qn \ 


duPn 

V o / 

Q c n ,P^ satisfy the condition 
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L p 
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\ 
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\ 



lim 

h^o+ 



Qn 



\ P n j 
Qr 



Qn 


V P n ) 

Q{x) 
P(x) 



( dl-Qn \ 



d^Pn 

V o J 



(2.8) 



(2.9) 



where Q(x),P(x) are the continuum limit of Q n (0-> P n (0- The eigenvalue at A = has algebraic 
multiplicity six. 



3 Evans Function and stability 

We rewrite the eigenvalue problem (L — A)u = as a system of difference equations 

Y n+1 = A(A, n)Y n 



(3.1) 



It is known that there exist solution sets Y^,Y^,Y^ and Y^,Y^,Y^ to equation (3.1) such that 



-> exponentially fast as n — > — oo for i = 1, 2, 3 and 
i = 4,5,6. 

The Evans function associated with (3.1) is given by 



exponentially fast as n — > oo for 



n—l 



£(A) = (Y-AY+)/;Qdet(A(A,j)) 

3=0 

satisfies the following properties [43]: E(\) is analytic in A, E(X) = iff equation (3.1) has a bounded 
solution and the order of the zero is equal to the algebraic multiplicity of the eigenvalue. 



The adjoint problem associated with equation (3.1) is 

Z n+1 = [A(A,n)- 1 ]*Z r , 

and its solution satisfies: 



Sij i,j = 1,2,3,4,5,6. 



Assuming that E{X) / 0, these solutions furthermore have the property 



(3.2) 



exponentially fast as 



oo for % = 1, 2, 3 and 



exponentially fast as n — > — oo for i = 4, 5, 6. 



For the particular problem of CDNLS, we obtain the linearized system (3.1) with 
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\ 
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where 



1q = a{(3 - 2Q n {Q n+1 + Q n _i) - 1), 7p = a(/3 - 2P n {P n+1 + P n _i) - 1). 



The solutions of (3.1) about A = are the following 
Yjr(0)=Y£+(0) = 
Y T 2 r(0)=Y^+(0) = 



%Q n , 0, %(Q n - Q n _i)//i, 0, % P n , 0, %(P n - P n _i)//i, 

T 



T 



T 



0,Q n ,0,(Q n -Q n _i)//i,0,0,0,0 

Y3.-(0)=Y3.+(0) = [0,0,0,0,P n ,0,(P n -P n _i)/h,0 
Furthermore, define the relevant adjoint solutions as 

d^Qn-i - Qn)/h, 0, dt.Q ni 0, %(P n -l - P n )//i, 0, 9 5 P„, 

T 



T 



y2 _ 



0, {Qn-l-Qn)/h,Q,Qn, 0,0, 0,0,0 



0,0,0,0,(P n _i -P n )//i,0,P n ,0 



T 



Finally for the problem Y n = A(0, n)Y define the solutions: 



u r 



Yi'-(0), ^ = Y^-(0), i4 = Y3."(0) 



and let and it„ satisfy the relations 



4-zi = o, i4-zi = i, ^-zi=o 
t4-z£ = i, i4-z2=o, i4-z* = i 

7 



(3.3) 



(3.4) 



(3.5) 



v*-Z 3 n = l, u 5 n -Z 3 n = 0, u 6 n -Zl = 0. 

The equation of variation with respect to e for the stable W s and W u manifolds is given by the 
non-homogeneous problem 

Y n+ i = A(A,n)Y n + g n 



where { g ra } is a uniformly bounded sequence: 

g n = -e2h{ql+p 2 n ) 



( hq n \ 

hp n 
V Pn ) 



The distance between stable and unstable manifolds can be calculated by 

d £ (W s -W u ) = M(^h)ut 

where M(£, /i) is the Melnikov sum 

CO 

M(£,h)= ]T Sn-Z„+1 
n=— oo 

After the substitution of homoclinic solutions, the Melnikov function takes the form 

oo 

M(£,h)=-2h ]T (Q 2 n + Pn){Qnd^Q n + Pnd^Pn). 

n=— oo 

Let us consider 



Qisinh2Vy 



Then, the homoclinic solutions (2.7) are given by 



«2 



a2sinh2VF 



Qn(0 = uiQ(0, Pn(0=a2Q(0, where Q(£) = sech(Wn + £)■ 
It is clear that the Melnikov sum can be rewritten as 

h 



Using the Poisson summation formula we obtain 

M(£, h) = a w C M sm{2^/W) + 0(e" 2 ^) 

where 



1-2. i ~2\2 

= K + a 2 ) , 

„ . ~ N h , 2ir .. 4 7r 4 7T 3 . 
C„(»')=V( 17 )(- w + - ip) e- 



(3.6) 



(3.7) 



(3.8) 



In order to calculate the Taylor series expansion of the Evans function about the eigenvalue A = 0, 
we use the fact that A = has algebraic multiplicity six. One must therefore calculate d®E(0) (see 
Appendix): 



where 



and 



dtE(0) = ^^BxB^a. (3.9) 

OO -i p oo 

B l= E l(Q C n^Qn + P^P n )= - _ (Q 2 (x)+P 2 (x))dx + 0(h), (3.10) 

n=— oo 00 

00 a 1 Z" 00 

B 2 = E -TQnd u Qn = -^d u Q 2 (x)dx (3.11) 
n=— 00 00 

1 /'OO 

B 3= E ~^Pnd UJ P n = --d w P 2 (x)dx + 0(h). (3.12) 



Due to the fact that CDNLS system is Hamiltonian, the eigenvalues for the linearized problem will 
satisfy the relationship that if A is an eigenvalue, then so is —A and ±A*. Furthermore, for the CDNLS 
[0, 0, Q n , P n ] T remains an eigenfunction, and, in turn, [d^Qn, d^Pn, 0, 0] T remains a generalized eigen- 
function; hence A = is an eigenvalue with multiplicity four in the perturbed problem. Consequently, 
we need to compute d £ dfE(0). One can write 



d £ diE(0) 



d £ (Ylr A Y^+) A d 2 {Yl- A Y£+) A d 2 (Y^ A Y T 3 t >+) A Y 1 * A Y 2 n + A 



(0). (3.13) 



As above and using 

d £ (Y 1 ^AY 1 ^) = d^M(C,h)ul 

we obtain 

d £ 8{E(0) = -±B 2 B 3 ^M(e, h). (3.14) 
Combining the above results one has now an expansion of the Evans function about A = as 

E{\) = B 2 B 3 -^r [ jU%M(£, h) + ^A 2 Bi ] , (3.15) 

where M(£,h) is given by (3.8). 

The following theorem can then be stated. 

Theorem 1 Consider the stability of the solitary waves associated with the CDNLS for h > sufficiently 
small. The associated linear operator has four eigenvalues at A = 0. Furthermore, there exist only six 
eigenvalues near A = 0. If £ = 0, then the wave is linearly stable, and the two additional eigenvalues 
are purely imaginary and are given to lowest order by 



V A ^i w 

whereas if £ = W/2, then the wave is linearly unstable, and the two additional eigenvalues are real and 
are given to lowest order by 

A± = ± \^ea w C M 



4BiVF 
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The predictions of Theorem 1 were numerically tested in Figs. 1 and 2 for the site-centered and 
inter-site centered modes numerically. The solutions were constructed for different values of e, using a 
fixed point iteration scheme on a 100-site lattice. The starting point used was the analytically available 
solutions [46] at the limit of e = (i.e., the integrable limit). Upon convergence to the exact numer- 
ical solution for finite e, numerical linear stability analysis was performed to obtain the spectrum of 
linearization around the solitary waves. We observed that as soon as one deviates from the integrable 
limit, for this nonlinearly coupled case, the effective translational invariance of the discrete model is 
"broken" resulting in the bifurcation of the relevant pair of eigenvalues from the origin of the spectral 
plane. This bifurcation will occur along the imaginary axis for the site-centered mode of Fig. 1, while 
the eigenvalues will exit as real in the inter-site centered case of Fig. 2, in accordance with Theorem 1. 
Furthermore, the dependence of the eigenvalues will be essentially proportional to \fi as is quantified 
in the right panel of the relevant figures, which is consonant with the prediction of Theorem 1 above, 
according to which A oc e 1 / 2 . On the contrary, the two phase invariances of the two components, corre- 
sponding to the respective "mass" conservation laws norm of each field are preserved by the nonlinear 
coupling. As a result, the relevant 2 pairs of eigenvalues at the spectral plane origin, will remain at 
A = 0. 

We note in passing that another implication of Theorem 1 concerns the near continuum limit behavior 
of the relevant eigenvalues, according to which: 

for small h, which is consistent with the earlier findings of [43] (see also references therein about exponen- 
tially small splittings of eigenvalues in the presence of discreteness). Hence, approaching the continuum 
limit (for fixed e, and h — > 0), the relevant translational eigenvalue pair would approach A 2 = 0, in 
direct analogy to its scalar case [43] (given this analogy, we do not examine this case further). In the 
line of this qualitative analogy (the quantitative details and relevant constants differ due to presence of 
two components here), Theorem 1 can be parallelized to Theorem 3.4 of [43]. 

A further note can be added regarding possible variations of the scattering lengths, which is rele- 
vant e.g., to our motivating example of two hyperfine states of 87 Kb [13, 14]. Slight deviations of the 
scatterings lengths from their unit values as well as possible changes of the coupling constants among 
the different components do not change the relevant eigenvalue count (since they do not affect the 
symmetries of the problem) or the over-arching stability conclusions for on-site versus inter-site modes. 

4 Coupled DNLS with Linear Coupling 

In this section, we consider a system of CDNLS equations with linear coupling: 

\u n = -A 2 u n - (\u n \ 2 + \v n \ 2 )(u n+1 + u n _i) + Sv n 

\v n = -A 2 v n - (\v n \ 2 + \u n \ 2 )(v n+ ± + v n -i) + Su n . (4.1) 
We are interested in stationary solutions of the system of equations (2.1) and make the ansatz 

u n {t)=q n ^\ v n (t)=p n e-^, (4.2) 
with a uniform rotation frequency to. Upon setting 

r n = (q n ~ q n -i)/h, s n = (p n - p„_i)//i, 
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Figure 1: The left panels of the figure show the branch of stable, site-centered solitons stemming from 
the continuation of the integrable case of e = 0: the top panel shows the norm of the solution as a 
function of e; the middle shows the discrete soliton spatial profile for e = 1 and the bottom panel shows 
the spectral plane of the linearization eigenvalues around the solution (again for e = 1). The right 
panels show the evolution of the eigenvalue bifurcating from the origin of the spectral plane (due to 
the breaking of the "effective" translational invariance of the integrable case. The trajectory of this 
imaginary eigenvalue pair is shown as a function of e in linear (top panel) and log-log (bottom panel) 
plot. The best fit power law is shown by solid line and has the exponent p = 0.53. 




1 e 




Figure 2: The same features are shown as in Fig. 1 but for the inter-site centered localized mode, where 
the relevant translational eigenvalue pair becomes real rendering the configuration unstable. The best 
fit exponent is p = 0.54 in the right panel showing the real eigenvalue pair as a function of e. 
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the steady-state problem for CDNLS can be formulated as 



where 



Q.n+1 = (ex - l)q n + r n h + 5h 2 p n x 

r n+ i = (c%-2)y +r n + 5hp n x 

Pn+i = (cx - l)p n + s n h + 5h 2 q n x 

s n+1 = (cx - 2) j- + s n + Shq n x 

1 



(4.3) 



c = 2 — toh and x = 



l + h?{q 2 n +piy 

The equations (4.3) are written so that, when 5 = 0, they revert to the steady-state problem for 
integrable CDNLS. 

The equation of variation with respect to e for the stable W s and W u manifolds is given (as in the 
CDNLS (2.1) by the non-homogeneous problem 

Y n+ i = A(A,n)Y n +g„, (4.4) 

where { g n } is a uniformly bounded sequence: 

/ hpnX \ 
gn = 6k hq nX 

V q n x- ) 

The distance between stable and unstable manifolds can be calculated by 

d s (W s - W u ) = M 2 (lo, h)u A n . (4.5) 

M.2(ijJj h) is the Melnikov sum 

oo 

M 2 (uJ,h) = §n • 

n=— oo 

where the factor x in g is eliminated because we examine what happens for small h and 

-i T 

du,{Qn-l ~ Qn)/h, d^Qn, ^(P n _l - P n )/h, d w P n 

the second "growth mode" is obtained by solving the inhomogeneous equation. 
It is clear that the Melnikov sum can be rewritten as 

oo ^oo 

M 2 (w, h) = h5d w p nQn = 5d w / P{x)Q{x)dx + 0(h), 

where Q n ,Pn are defined in (2.7). Using the continuum limit of P and Q, one can approximate M 2 
close to the continuum limit as: 

a\OL 2 



M 2 = - 



a\ + a 2 



<5(-^)- 1 / 2 . 



(4.6) 



Combining the above results one has now an expansion of Evans function about A = 0, similarly to 
Eq. (3.15), involving the derivatives dsd x (0) and d®E(0), as the leading terms; the details are left to 
the interested reader, being quite similar to those of the previous section. One then has the following 
theorem: 
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Figure 3: The left panels are similar to those of Fig. 1, but are now for the site-centered mode in 
the case of linear coupling between the components. The middle and bottom left panels show the field 
profile and the stability for the case of 5 = 1. The right panel shows the eigenvalue bifurcating from the 
origin as soon as 5 becomes non-zero. The circles denote the full numerical results, while the solid line 
denotes the curve Aj = 28, which approximates very well the numerical result. 



Theorem 2 Consider the stability of the solitary waves associated with the CDNLS (4-1) for h > 
sufficiently small. The associated linear operator has four eigenvalues at A = 0. Furthermore, there 
exists an additional pair of eigenvalues near A = which are purely imaginary and directly proportional 
to S. 

We have tested the above predictions numerically and have found them to be in good agreement with 
the computational results. In particular, the site-centered and inter-site centered modes are respectively 
shown in Figs. 3 and 4 for the case of linear coupling (in order to compare/contrast them with those of 
nonlinear coupling in Figs. 1-2). The main thing to notice in this case is that the eigenvalue bifurcation 
is not only linear in its dependence of the linear coupling parameter a, but is furthermore along the 
imaginary axis for both the site-centered and inter-site centered solutions. The latter implies the absence 
of instability in such linearly coupled cases. It is important to highlight here a key difference between 
the linear and nonlinear coupling case, to which this difference in stability properties can be attributed. 
In the latter case, examined previously, the effective translational invariance of the integrable limit was 
"destroyed" upon the action of the nonlinear coupling, while the relevant phase invariances remained 
intact. On the contrary, in the linearly coupled case, the translational invariance remains present, while 
one of the phase invariances is destroyed, as it is now only the sum of the squared I 2 norms that is 
conserved (rather than each individual one of them). As a result, the source of the bifurcation, and 
hence the ensuing stability behavior is different in the two cases. This is also implicitly mirrored in 
the similar stability behavior of the site-centered and inter-site centered modes. Finally, let us make, 
in passing, another important observation for the linear case. For the setting examined above with 
solutions identical between the two components i.e., with q n = p n , it is straightforward to note that 
the exact solution is analytically available for all values of the linear coupling 5, as the latter merely 
renormalizes W, by replacing uj with w — 5 in the relevant expression. Furthermore, an alternative 
way to see that the case with the linear coupling should lead to linear dependence of the near A = 
eigenvalue pair as a function of 5 is the following. Consider the linear transformation [49] 

u n \ _ I cos(5t) —isin(St) \ / u n \ 
v n J y -ism(5t) cos(5t) J \ v n J 
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Figure 4: The same features are shown as in Fig. 3 but for the case of the inter-site centered mode. The 
middle and bottom left panels are for 5 = 0.5. The right panel shows the eigenvalue bifurcating from 
the origin. The circles denote the full numerical linear stability results, while the solid line denotes the 
curve Aj = 25. 

It is interesting to observe that the equations satisfied by u n and v n are those of the original model i.e., 
without the linear coupling. Hence, the linear coupling can be "factored out" through this transforma- 
tion, which is motivated by first-order vector systems of ordinary differential equations (and respected 
by our Manakov-type nonlinearity) . The eigenvalues of the transformation matrix are exp(i5t) and 
exp(— i5t), which in turn suggests a linear dependence of the bifurcating pair of eigenvalues on 5 in 
agreement with our numerical results of Figs. 3-4. 

5 Conclusions 

In conclusion, in this paper we have developed the Evans function methodology for discrete nonlinear 
Schrodinger equations in the vector case. We have used as our starting point the integrable discretization 
of [44] and have introduced nonlinear, as well as linear perturbations to it breaking different kinds of 
invariances. As a result, pairs of eigenvalues, corresponding respectively to these invariances (of the 
linearization around solitary wave solutions of the equation) have moved away from the origin of the 
spectral plane. These pairs have been analytically tracked via the zeros of the Evans function and 
have been found to yield critical information about the stability of the solutions in the non-integrable, 
perturbed case. 

In the nonlinearly coupled case, an eigenvalue pair corresponding to translations of the solitary waves 
has been found to bifurcate from 0, leading to either a stable (site-centered) or unstable (inter-site cen- 
tered) solution. On the other hand, in the linearly coupled case, a pair corresponding to the phase 
invariance of the waves has been found to bifurcate linearly from the origin. In both settings, the com- 
puted results based on numerical existence and linear stability methods have successfully corroborated 
the analytical predictions. 

It would be interesting to extend the methodology presented herein to the, very intensely studied 
in recent years, context of photorefractive materials [50, 51]. In the latter case, the nonlinearity is of 
the saturable type (coinciding with the ones studied above for low intensities, but having very different 
-quasi- linear- behavior at high intensities). Such studies and the corresponding extensions to higher- 
dimensional settings are currently in progress and will be reported in future publications. 
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A Calculation of d 6 x E{0). 

Here, we calculate the d®E(0). Following Kapitula [47], [43] one has 



3$E(0) = 120 



dl{Y^~ A Y£+) A dl(Y*r A Y T 2 t >+) A ^(Y*" A Y£+) A Y};+ A Y 2 '+ A Y*" 



For i = 1, 2, 3 we have 



dlY^O) = A(0,n)<9 2 Y^(0) + 23 A A(0,n)3 A Yj; ± (0), 



(0). (A.l) 



(A.2) 



where 



d X Y^(0) = [0,Q c n ,0,(Qn-Qn-l)/h,0,Pn,^(Pn-Pn-l)/h 

d x Y^(0) = [d u Q n ,0,d u (Qn-Q n -i)/h,0,0,0,0,o] T 
dxY^iO) = [0,0,0,0,^P n ,0,^(P n -P n _i)//i,0, ] T 



(A.3) 



Set gi = d A A(0,n)d A Y^(0) for i = 1,2,3. Upon using the variation of parameters to solve equation 
(A.2), one has that 



flS(Y£-AY£+)(0) 
M"A^ + )(0) 
^(Y«.-AY 3 fl .+)(0) 



2 



2 / 

^(«n S Sn • Z n+1 + Cl«J, + C 2 U n + C 3 *4 
n=— oo 

2 / 00 \ 

Z • K+l + + C 5 U n + C 6 U 3 n J 

^ n=— oo ' 

(«b E Sn ■ Z n+1 +C7< + C 8 ^ + c 9 n^y 



(A.4) 



Here ci, ■ • • , eg are constants. The above equations can be simplified as follows 



flSOft-AYi'+XO) 



oo 

7 ") I 2 f ( Qn^Qn + ^^P„ ) + C X u\ + Cfcli* + 

?i=— oo 
oo 



U, 



C 3 U r 



^(YfAYf)(0) = 
OT-A^+)(0) = 



2 
2 



U ™ Z ~ ^ QnduQn + CAU l n + C^u\ + C 6 ?4 J 
n=— oo ' 

(°° Q! \ 



(A.5) 
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